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A general continuum theory is developed for ion intercalation dynamics in a single crystal of a 
rechargeable battery cathode. It is based on an existing phase-field formulation of the bulk free 
energy and incorporates two crucial effects: (i) anisotropic ionic mobility in the crystal and (ii) 
surface reactions governing the flux of ions across the electrode/electrolyte interface, depending 
on the local free energy difference. Although the phase boundary can form a classical diffusive 
"shrinking core" when the dynamics is bulk-transport-limited, the theory also predicts a new regime 
of surface-reaction-limited (SRL) dynamics, where the phase boundary extends from surface to 
surface along planes of fast ionic diffusion, consistent with recent experiments on LiFeP04. In the 
SRL regime, the theory produces a fundamentally new equation for phase transformation dynamics, 
which admits traveling-wave solutions. Rather than forming a shrinking core of untransformed 
material, the phase boundary advances by filling (or emptying) successive channels of fast diffusion 
in the crystal. By considering the random nucleation of SRL phase-transformation waves, the theory 
predicts a very different picture of charge/discharge dynamics from the classical diffusion-limited 
model, which could affect the interpretation of experimental data for LiFeP04. 



I. INTRODUCTION 

LiFeP04 is widely considered to be a promising cath- 
ode material for Li-ion rechargeable batteries. The high 
practical capacity and reasonable operating voltage of 
the material, along with its nontoxicity and potential low 
cost, make it well-suited for large-scale battery applica- 
tions P, 0, [1] . Unlike many other cathode materials that 
increase their Li concentration in a continuous solid so- 
lution, Li2;FeP04 only exists for a; « and x sa 1 Q 
and charges or discharges Li by changing the fraction of 
phase with x « and a; « 1. 

The transport and phase separation properties of 
LiFeP04 have been studied extensively by atomistic sim- 
ulations 0, @, 0j Q- First-principles calculations have 
shown that Li diffusion in the bulk FeP04 crystal is 
highly anisotropic d, 0] • Li is essentially constrained to 
ID channels in the (010) direction arranged in layers that 
form the crystal, as depicted in Fig. [1] The lattice mis- 
match at the FeP04/LiFcP04 phase boundary is signifi- 
cant (5% in the x direction, shown in Fig. [T|), and recent 
work has investigated the differences in clastic properties 
between the lithiated and unlithiated material [9| . Atom- 
istic simulations have also suggested that electrons in the 
crystal may diffuse as small, localized polarons confined 
to planes parallel to the Li channels [l3, El, 113 ■ 

Recent experiments have confirmed the anisotropic 
transport and phase separation of Li in single crystal 
LiFeP04 [H, [ll, [iBl- Moreover, detailed microscopy 
in these studies has revealed that the FeP04/LiFeP04 
phase boundary is a well-defined interface that extends 




(010) 



FIG. 1: Schematic of plate-like single crystals of LiFeP04. Li 
is confined to ID channels in the y direction, and channels 
are stacked in layers parallel to the yz plane, indicated by 
the dotted lines. Typical dimensions of single crystals are 
2 X 0.2 X 4 ^m in the x, y, z dimensions, respectively [l3l |. 
For each direction, the corresponding space group Pnma axis 
and Miller index are shown in parentheses. 



through the bulk crystal to the surface. In experiments, 
the phase boundary has a characteristic width of several 
nanometers on the surface although this width is 
probably broadened by experimental resolution, and Li 
insertion and extraction seem to be concentrated in this 
region, with negligible transfer occurring in cither the 
FCPO4 or LiFcP04 phases. Notably, the phase boundary 
moves orthogonally to the direction of the surface flux, in- 
dicating that as Li insertion (extraction) proceeds, layers 
of the ID channels are progressively filled (emptied). The 
observance of surface cracks and their alignment with 
the phase boundary [l^ also reinforces the view that 
the FeP04/LiFeP04 lattice mismatch plays an impor- 
tant role in the electrochemical function of the material, 
as may the associated stress field [l6j . 
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In light of the understanding gained from these atom- 
istic and experimental studies, the continuum theory of 
transport and phase separation in LiFeP04 merits re- 
newed attention. The prevailing "shrinking core" model 
can be traced to a qualitative picture accompanying the 
first experimental demonstration of the material as an 
intercalation electrode This model assumes a grow- 
ing shell of one phase surrounding a shrinking core of 
the other phase, with the shell and core phases deter- 
mined by the direction of the net Li flux: a LiFeP04 
shell surrounds an FeP04 core during Li insertion (bat- 
tery discharging); an FeP04 shell surrounds a LiFeP04 
core during Li extraction (battery charging). It is im- 
portant to note that the boundary between the shell and 
core phases is entirely contained within the bulk of the 
material and moves parallel to the direction of the Li flux, 
in contrast to the observations of the experiments cited 
above. 

The current state of mathematical modeling of ion in- 
tercalation is based on the shrinking-core concept with 
some further simplifying assumptions. In earlier work, a 
simplified version of the model was mathematically for- 
mulated by Srinivasan and Newman p7| and incorpo- 
rated into an existing theory for transport in composite 
cathodes [l^. In their model, FeP04 is treated as a 
continuous, isotropic material, and Li is inserted and ex- 
tracted uniformly over the surface of a spherical FCPO4 
particle. The phase boundary is defined as where the 
compositions LieFeP04/Lii_eFeP04 coexist, with e <C 1 
specifying the equilibrium composition between the Li- 
poor and Li-rich phases, and no nuclcation constraints 
are included. Only Li transport in the shell is considered 
and modeled by an isotropic, constant diffusivity diffu- 
sion equation, while the velocity of the phase boundary is 
prescribed by a mass balance across the boundary. Thus, 
for Li insertion (extraction), diffusion in the growing shell 
occurs between the surface Li concentration and 1 — e (e) . 
The value of e is set as a parameter in the numerical so- 
lution of the model. 

In this paper, we propose a more general contin- 
uum theory for ionic transport and phase separation 
in single-crystal rechargeable battery materials, focus- 
ing on the special case of LiFcP04. Our theory ac- 
counts for anisotropic Li diffusion in the bulk as well 
as the formation and dynamics of the FeP04/LiFcP04 
phase boundary, driven by surface reactions at the elec- 
trolyte/electrode interface. We utilize an existing phase- 
field model for the free energy of the system to calculate 
the Li chemical potential [l9[; the bulk transport equa- 
tion and surface reaction rates for Li are then derived 
in terms of this chemical potential. The phase-field ap- 
proach provides a sound thermodynamic basis for study- 
ing the system and also directly connects our theory to 
atomistic modeling, since first-principles simulations can 
accurately compute the Li chemical potential [l^, [2l| . 

We first develop a general model that encapsulates var- 
ious transport regimes in different limits of the character- 
istic timescales for bulk diffusion and surface reactions. 



as presented in Fig. [2l The characteristic timescales for 
bulk diffusion to and surface reactions are 



to 
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where L is the lengthscale over which diffusion occurs, D 
is the diffusivity, and k is the surface reaction rate. The 
dimensionless ratio of these timescales is the Damkohler 
number 



Da 



tR 



IT' 



(2) 



Therefore, an isotropic bulk transport limited (BTL) pro- 
cess, where bulk diffusion in all directions is much slower 
than surface reactions, is characterized by 



Da > 1. 



(3) 



In this regime, the phase boundary is entirely contained 
within the material and moves along the direction of the 
Li flux, as shown in Fig. [5^. Anisotropic BTL phase 
transformation in LiFcP04 is depicted in Fig. [2]d. Bulk 
diffusion in the x and z directions is negligible, and bulk 
diffusion in the y direction is much slower than surface 
reactions. Consequently, 



Dax, Daz > Day > 1, 



(4) 



where Da^. = kL^ / D^ with D^ the diffusivity in the x di- 
rection, and similarly for Da.y and Da^- In this regime, 
the phase boundary is still contained within the bulk, 
but Li is confined to ID channels in the y direction. The 
anisotropic BTL model is a generalization of the shrink- 
ing core model developed by Srinivasan and Newman 
17|. However, motivated by the experiments and simu- 
lations described above, we focus on a different transport 
regime where surface reactions are much slower than dif- 
fusion in the y direction but much faster than diffusion 
in the x and z directions. Thus, we study the limit 



Dax, Daz '> !:$> Da 



(5) 



illustrated in Fig.^:. The phase boundary in this regime 
extends through the bulk of the material to the xz sur- 
faces. Moreover, the surface reaction rates of Li trans- 
fer are concentrated at the phase boundary such that 
each channel is almost completely lithiatcd (delithiated) 
before insertion (extraction) progresses to an adjacent 
channel. In such an anisotropic surface reaction limited 
(SRL) process, our general model reduces to a fundamen- 
tally new equation governing the transport and phase 
separation dynamics. Analysis and numerical solutions 
of this equation show that it qualitatively reproduces the 
features of the LiFeP04 system observed in experiments. 

In this initial effort, we neglect the possibility of charge 
separation and assume that electrons are freely available 
in the material to compensate the charge of Li"'', although 
the presence of multiple diffusing and migrating species 
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FIG. 2: Transport models obtained in different limits of the characteristic timescales for bulk diffusion and surface reactions. 
Figures show xy cross sections of spherical and plate-like single crystals during Li insertion, after phase nucleation has occurred. 
Lithiated portions of the crystal are shaded, and points outside particles represent flux of Li ions across the electrode/electrolyte 
interface (shown only for spherical particles). The FeP04/LiFeP04 phase boundary is denoted by the dashed line, and arrows 
indicate movement of the boundary as insertion proceeds, (a) Isotropic bulk transport limited, (b) Anisotropic bulk transport 
limited, (c) Anisotropic surface reaction limited. 



(interacting through an electrostatic potential) can be 
modeled as an extension of the general framework we 
present. We also avoid an explicit treatment of stress 
around the FeP04/LiFeP04 interface; instead, we note 
that a term in the phase field formulation may serve as an 
approximation of the energy due to the lattice mismatch. 



II. GENERAL MODEL 

A. Phase field formulation 

We follow the conventional Cahn-Hilliard formulation 
[2^ that has been previously applied as a phase field 
model for bulk transport in LiFcP04 [l^. The total free 
energy of the system is expressed as a functional of the 
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local Li concentration 



I [/(c) + (/^/2)|Vc|2] dr, 



(6) 



where c is the dimensionless, normalized Li concentra- 
tion (0 < c < 1), /(c) is the homogeneous free energy 
density, and K is the gradient energy coefficient that 
represents the energy penalty for maintaining concen- 
tration gradients in the system. The lattice mismatch 
at the FeP04/LiFeP04 phase boimdary coincides with 
the concentration gradient, so the gradient penalty in ([6|) 
can be regarded as approximating related contributions 
the free energy. Phase-field models have also been devel- 
oped, which include the long-range elastic contributions 
to the free energy [H, , but here we restrict ourselves 
to the simpler formulation above, to allow us to focus on 
other effects, namely anisotropic transport and surface 
reactions. In this spirit, we also ignore any surface con- 
tributions to the free energy, such as the tension of the 
electrode/electrolyte interface. 

The homogeneous, bulk free energy density takes the 
form of a regular solution model 

/(c) = ac(l -c)+ pfcr[clnc+ (1 -c)ln(l - c)] , (7) 

where a is the average energy density (in a mean field 
sense) of the interaction between Li ions, p is the number 
of intercalation sites per unit volume, k is Boltzmann's 
constant, and T is the temperature. The first term in ([7]), 
the enthalpic contribution to the free energy, promotes 
separation of the system to c = or c = 1, while the sec- 
ond term, the entropic contribution, favors mixing of the 
system. Therefore, the strength of the phase separation 
is characterized by the dimensionless ratio pkT/a. The 
chemical potential of Li in the FeP04 host is calculated 
as the variational derivative 



5c' 



(8) 

= p - KV^c, (9) 
where we define p. as the homogeneous chemical potential 



dl 
9c' 

a(l - 2c) + pkT In 



1 -c 



(10) 
(11) 



While in general, the two phase compositions in equi- 
librium across the miscibility gap are determined by the 
common tangent construction, in the symmetric free en- 
ergy / of ([7]) these compositions correspond to p ~ 0. 
These roots cannot be found analytically from pT|) . 
but asymptotic approximations in the small parameter 
pkT/a can be obtained; two term expansions are 



e pi'T [1 + 



pkT 



_g pkT 



c+~l 



(12) 
(13) 



where c_ is the root near c = and c+ is the root near 
c = 1. As may be expected, c± approach the concen- 
tration extremes exponentially in a/{pkT). While nucle- 
ation may be required to form a second phase for compo- 
sitions in the miscibility gap, spontaneous phase separa- 
tion occurs when the composition is within the spinodals. 
The spinodals correspond to the zeros of the curvature 
of the free energy and can be determined from ([7]) as 



1 ± ^1 - 2pkT/a 



(14) 



We observe that a > 2pkT is required for distinct, phys- 
ically meaningful spinodal compositions. 



B. Anisotropic bulk transport 

As described above, Li migration in the bulk crystal is 
confined to ID channels in the (010) direction, labeled as 
y in Fig. [2l Some diffusion in other directions may oc- 
cur due to defects in the crystal lattice or cracks caused 
by the FeP04/LiFeP04 lattice mismatch, as have been 
observed experimentally [l^, [l^. Experiments [H, [l3| 
have found that layers of stacked ID channels in the z 
direction are progressively filled or emptied as the phase 
boundary moves across the layers in the x direction, in- 
dicating that transport in the z direction is faster than 
transport in the x direction. The anisotropic Li flux can 
be written as 



-cBV/i, 



where 



'11 






















^33 



(15) 



(16) 



is the mobility tensor for an orthorhombic crystal, and 
the indices 1, 2, 3 correspond to the a;, y, z directions, 
respectively. The diffusivity tensor is determined from 
the mobility by the Einstein relation D = fcTB. Based 
on the discussion above, we expect that bu ^ 633 ^ 622- 
We also note that D can, in general, be position or con- 
centration dependent. In fact, for some other intercala- 
tion materials that do not phase separate, first principles 
calculations have shown that the diffusivity is a highly 
nonlinear function of the Li concentration, varying by 
several orders of magnitude over the composition range 

With the ionic fluxes thus defined, the dynamics of the 
concentration profile is governed by the conservation law 



dc ^ . ^ 
P^+V.j^O, 



(17) 



where the factor of p is needed for dimensional consis- 
tency. 
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Surface reactions 



The net rate of Li insertion is therefore 



We assume that Arrhenius kinetics govern the insertion 
and extraction rates of Li across the electrode/electrolyte 
interface. The activation energies for these interfacial re- 
actions are defined as the change in the chemical poten- 
tial of Li across the interface. We further assume that 
the chemical potential of Li in the FeP04 host, fi given 
by ©J is valid everywhere on the electrode surface where 
Li transfer occurs, with the appropriate modification for 
the Laplacian term. By also using the bulk chemical po- 
tential for ions at the crystal surface, we are neglecting 
the possibility of any variation in the chemical potential 
at the electrode/electrolyte interface, e.g. due to surface 
orientation or surface curvature, as noted above [1^. 

With these assumptions, the local Li insertion rate is 



,(Ale-M)/(pfcT) 



(18) 



where /cins is the insertion rate coefficient, and Cg and fie 
are the concentration and chemical potential of Li in the 
electrolyte, respectively. Note that since Ce is expressed 
as a dimensionless filling fraction, i?ins and fcjns have di- 
mensions of inverse time. 

In a more complete battery model, Cg and fi^ would be 
determined by solving the appropriate transport equa- 
tions for Li in the electrolyte. However, as we focus here 
on Li transport in the electrode, we ignore variations in 
the electrolyte and take Ce and fie as constants. Our for- 
mulation therefore describes potentiostatic, or constant 
chemical equilibrium, conditions in the electrolyte, with 
the interfacial transfer of Li as the rate limiting process. 
Absorbing Cg into fcjns gives 



aKV^c 



(19) 
(20) 



where a = 1/ [pkT) , and i?ins is the homogeneous inser- 
tion rate 



-"'ins '^ms 



1 - C 



Similarly, the extraction rate is 



cxt "-cxt 



where i?ext is the homogeneous extraction rate 

^2 



R 



cxt "'CXt 



1 - c 



„Q[a(l-2c)-A<e] 



(21) 



(22) 
(23) 



(24) 



In contrast to i?insi the proportionality of i?cxt on c can- 
not be neglected and breaks the symmetry about c = 1/2. 



R R'ms Rc->ct 7 



1 -C 



1 - C 



^a[tJ.,-a(l-2c}+KV^c] 



„a[a(l-2c)-^e-_fS'V^c] 



(25) 
(26) 

(27) 



The boundary conditions for (jl7p on the crystal surface 
express mass conservation at the electrode/electrolyte in- 
terface. 



n • J 



-PsR, 



(28) 



where n is the unit normal vector directed out of the crys- 
tal, and ps is the number of intercalation sites per unit 
area, which depends on the orientation of the surface. 
Consistent with our neglect of surface excess chemical 
potential, we neglect the possibility of a surface flux den- 
sity js, whose surface divergence Vs • js would appear as 
an extra term on the right hand side of 



III. 



SURFACE-REACTION-LIMITED PHASE 
TRANSFORMATION 

A. Depth-integrated dynamics 



We now develop a special limit of the general model 
that describes SRL phase-transformation in LiFeP04. 
We assume that fast diffusion in the y oriented ID chan- 
nels rapidly equilibrates the bulk Li concentration to the 
surface concentration. A depth-averaged concentration c 
on the xz surface can therefore be defined as 



1 



Ly [X, Z) 



c{x,y,z,t)dy, 



(29) 



where Ly{x, z) is the depth of the crystal in the y direc- 
tion, from surface to surface. 

By depth- averaging the bulk transport equation (jl7p . 
using the boundary condition (|28p . we find that the dy- 
namics of c are governed solely by the surface reaction 
rate i?, acting as a source term on the xz surface. 



pLy{x, z) \ dc 



2ps{x, z) J dt 



R{c,V'c), 



(30) 



where Ps{x, z) is the number of atoms per unit area, de- 
pendent on the local orientation of the crystal surface. 
As noted above, a surface diffusion term could also be 
explicitly added to ([50]) . but we shall see that the reac- 
tion rate R already produces weak xz diffusion due the 
gradient penalty K, which suffices to propagate the phase 
boundary along the surface. 

Equation (j30p describes a fundamentally different type 
of phase transformation dynamics. From a mathematical 
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point of view, the striking feature is that the Laplacian 
term V^c appears in a nonlinear source term, as opposed 
to the additive quasiUnear term in classical reaction- 
diffusion equations. We are not aware of any prior study 
of this type of equation, so it begs mathematical analysis 
to characterize its solutions. 

Here, we begin this task by making some simplifying 
assumptions, which allow us to highlight new nonlinear 
wave phenomena described by pO)) . As noted above, ex- 
periments indicate that layers of ID channels along the 
yz plane are progressively filled or emptied as Li transfer 
proceeds, so it is natural to neglect concentration vari- 
ations in the z direction as a first approximation, for a 
planar phase boundary spanning the crystal. We also ne- 
glect depth variations, Ly = constant, and assume con- 
stant surface orientation, ps = constant, which corre- 
sponds to the common case of a plate-like crystal. 



B. Dimensionless formulation 

A dimensionless form of the model, suitable for math- 
ematical analysis, is found by scaling each variable to its 
natural units. The Li concentration is already expressed 
as a dimensionless filling fraction per channel c, which 
will depend on the dimensionless position x = x/ L and 
time t = t/r. Position is scaled to a length L, which 
characterizes the size of the crystal surface along which 
the SRL phase transformation propagates. The natural 
time scale from (PU]) is 



pLy 



2psku 



(31) 



which is the time required for the insertion reaction to fill 
a single fast-diffusion channel in the crystal (from both 
sides). Note that this time scale is proportional to the 
depth of the crystal, Ly. 

There are four dimensionless groups which govern the 
solution. The first is the ratio of reaction-rate constants 



A is an atomic length scale (1 A- 10 nm) much smaller 
than the crystal size (10 nm - 10 /im), the parameter A 
is typically small and lies in the range 10^"^ < A < 1. 

With these scalings, the SRL phase-transformation 
equation (j30[) takes the dimensionless form 



dc 
dt 



1 - c 



1 - c 



gAo-a(i-2e)+A''|^ 



^a(l-2E)-/i,-A^ 



(35) 



We will study solutions to this new nonlinear partial dif- 
ferential equation in the following sections, but we al- 
ready can gain some insight by considering the limit of 
a sharp phase boundary, A ^ 1, as discussed above. Ex- 
panding (j35p for small A, we obtain a reaction-diffusion 
equation at leading order. 



dc 

at 



dx^ 



{Rins — Rcxt), (36) 



where R-ms and i?oxt are the dimensionless homogeneous 
reaction rates 



R'ms — 



1 



1 - C 



„Ac-a(i-2E) 



„a(l-2e)-Ae 



(37) 
(38) 



Thus, we see that ([55)) has a direct analog to a reaction- 
diffusion equation with a weak, concentration dependent 
diffusivity and nonlinear source term, in the appropri- 
ate physical limit of an atomically sharp phase bound- 
ary. The detailed structure and dynamics of the phase 
boundary, however, must be obtained by solving the full 
equation ([35]) . Representative plots of the homogeneous 
net insertion rate i?i,is — Roxt are shown in Fig. [31 



C. Traveling waves 



(32) 



which measures asymmetry in the extraction and inser- 
tion reaction kinetics. By scaling energy density to the 
thermal energy density pkT, we arrive at three more di- 
mensionless parameters: 



pkT' 



Me 



^J■e 



pkT' 



and 



K 



pkTL^ 



(33) 



(34) 



The latter formula makes it clear that the natural length 
scale for the phase boundary thickness, set by the gradi- 
ent penalty in the free energy, is A = a/ K/{pkT). Since 



We seek traveling wave solutions of ([55]) . which physi- 
cally correspond to waves of phase transformation, prop- 
agating through the FeP04 crystal with steadily trans- 
lating concentration profile. Since c is a depth averaged 
concentration, these composition waves extend through 
the bulk material to the surface, and they move parallel 
to the surface. Substituting the traveling-wave ansatz 



c(i, t) = g{x — vt), 



(39) 



into ()35p . where v is the unknown constant velocity, and 
rewriting the resulting ordinary differential equation as 
a first-order system gives 



9 
h' 



A2 



-vh + ^{vhy + Aug 



(40) 
(41) 
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FIG. 3: Dimensionless homogeneous net insertion rate, (a) 
/ie =0, ^ = 1, 5 = 1, . . . , 5 as labeled in figure, (b) a = 5, 
ft = 1, /ie = —2, . . . , 2 as labeled in figure. 



where g and ft. are functions of C = i - 
points of the system are given by 



a(l - 2g) + hi 



1-.9 



wi. The stationary 

ft 0, (42) 
lnr = 0, (43) 



where we define another dimensionless constant F 



The solutions of (j43p correspond to the roots of the 
spatially homogeneous source i?ins — ^cxt and therefore 
determine the equilibrium phase compositions of the sys- 
tem. Explicit expressions for the solutions are not possi- 
ble, though at least one solution must exist, since the left 
hand side changes sign over the interval < g < 1 . How- 
ever, phase separation into Li-poor and Li-rich phases 
requires three solutions 51 < (?2 < .93, where g\ and gz 
are the equilibrium compositions bounding the moving 
wavefront, and g^ is an unstable intermediate composi- 
tion. We observe that the phase compositions are in- 



dependent of the dimensionless gradient penalty A^. At 
the threshold where transitions from one solution to 
three, the solutions and extrema of the equation coincide. 
The extrema occur at the compositions 



9± = 



(1 + 45) ± Viea^ - 40a-|- 1 



(44) 



and as two distinct extrema are needed for three solutions 
of (|43p . phase separation requires d > 5/4 -f ■\/3/2 w 
2.47. If the minimum d is exceeded, P5|) can be used 
to compute the rate coefficients and electrolyte chemical 
potential, expressed through the combination F, that will 
make either 17-1- the critical composition at threshold. For 
strongly phase separating systems, such as LiFcP04, we 
find asymptotic approximations of the solutions of (j43p 
in the small parameter 1/5. Two term expansions for 
each root are 



(45) 




52 



53 



1 

2 

1 ^ 





'1 5 " 




5 ^ 252 







(46) 
(47) 



The existence of traveling waves and the selection of 
the wave velocity in a phase separated system can be un- 
derstood by a linear stability analysis of ((40|) - (|4T|) about 
the three stationary points (g^, 0) for i = 1, 2, 3. We find 
that (51, 0) and (53, 0) are saddle points for all velocities, 
and (32,0) is either a stable node or stable spiral, de- 
pending on the velocity. Monotonic wavefronts between 
the equilibrium Li-poor and Li-rich phases correspond to 
trajectories in the (g, ft) phase space that connect (51, 0) 
and (53,0), bypassing {g-i^ 0). Following the continuity 
arguments presented in [27|, there is a unique velocity 
V such that the orientation of the eigenvectors at (51 , 0) 
and ((73, 0) allow a single trajectory joining these points. 
Thus, for a given set of parameters, all fully developed 
waves in a system propagate at the same velocity. 

A rigorous mathematical analysis of the traveling wave 
solutions of ([50]) . including their formal existence, stabil- 
ity, and velocity is beyond the scope of this work. Al- 
though analytical methods for studying traveling waves 
in parabolic systems are available |28| . they are usually 
developed for systems where there is a diffusion term plus 
a source independent of derivatives of the solution, as in 
To the best of our knowledge, (|30p represents a 
different type of equation admitting traveling wave so- 
lutions, where the curvature dependence of the source 
precludes the need for an explicit diffusion term. 

The nanoscale dimensions of the physical domain also 
complicate the analysis of ([50)1 . In other reaction dif- 
fusion equations, such as the Fisher equation, the wave 
velocity is determined by assuming an exponential de- 
cay of the leading edge of the wavefront as x ^ 00 [131 . 
A finite cutoff in the leading edge is known to signifi- 
cantly alter the velocity [2§| . Such cutoffs are present in 
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nanoparticles of LiFeP04, as the xz surface is bounded 
on the scale of the wave width. Moreover, the nanometer 
wave width describes the Li concentration across only a 
few atomic layers of the crystal, with each ID channel in 
the layer holding a single file of Li atoms. Therefore, c 
may be discontinuous for small particles. 



D. Wave propagation and nucleation 

We investigate the traveling wave solutions of the SRL 
phase transformation model by numerically solving ([35]) 
with an explicit finite difference method, second order 
in space and time. Representative phase transformation 
waves during Li insertion and extraction are shown in 
Fig. [H For these simulations, a = 5 and A = 1, similar 
to the values given in [l^. The choice of A corresponds 
to the dimensional lengthscales A = L = 10 nm. Al- 
though there is no available experimental or simulational 
data on the rate coefficients kms and k^xt , we expect that 
the insertion and extraction reactions occur on the same 
timescale and therefore assume k = 1. Insertion or ex- 
traction is forced by raising or lowering the electrolyte 
chemical potential to promote transfer in one direc- 
tion and inhibit transfer in the other; fl^ = 0.5 for the 
insertion process in Fig.[4K, and fie — —1 for the extrac- 
tion process in Fig. [^Js. The initial conditions for the 
insertion and extraction waves, denoted by dashed lines 
in the plots, are Gaussian fluctuations of the composi- 
tion representing nucleations of the lithiated and unlithi- 
ated phases, respectively. For the insertion in Fig. [3^, 
c{x,0) = 0.1 4- 0.8 exp(— a;^), and for the extraction in 
Fig. [Ij, c{x, 0) = 0.9 - 0.8 exp(-i2). 

The initial composition fluctuation rapidly develops 
into two wavefronts bounded by the equilibrium Li-poor 
and Li-rich phases gi and 33, respectively. The devel- 
opment of a fully formed wavefront involves both in- 
sertion and extraction. During the insertion process in 
Fig. [Hi; the maximum concentration of the initial fluctu- 
ation grows to (73, while the low concentration baseline 
decays to gi. Conversely, for the extraction process in 
Fig. [lb; the minimum concentration decays to gi, and 
the high concentration baseline grows to (73. Once fully 
developed wavefronts form, they propagate to the right 
and left with a constant velocity. We have verified that 
the velocity of a fully developed wavefront is constant for 
all times in the numerical simulations, as expected from 
the phase plane analysis described earlier. 

The dimensionless width w and velocity w of a fully de- 
veloped wavefront are determined by the parameters of 
the system. The main parameter controlling the width 
of the wave is A, since it contains the gradient energy 
coefficient K. As A is decreased, the energetic penalty 
for forming gradients in the concentration is lowered, re- 
sulting in sharper wavefronts spanning the equilibrium 
phase compositions. Additionally, we find that sharper 
wavefronts move at a slower velocity. The numerically 
computed dependence of the width and velocity on A is 
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FIG. 4: Numerically computed phase transformation waves 
during insertion and extraction. Dashed lines show initial 
conditions; solid lines show concentration profiles at uni- 
formly spaced times, with arrows indicating direction of wave 
propagation. For both insertion and extraction: a = 5, 
K = 1, A = 1. (a) Insertion wave with fie = 0.5, c(i, 0) = 
0.1 + 0.8exp(— i^). (b) Extraction wave with fie = — 1, 
c(f,0) = 0.9 - 0.8exp(-£^). 



shown in Fig. [5l It is apparent that both li and v scale 
linearly with A in the range of physical relevance. In fact, 
this linear scaling can be determined analytically by a di- 
mensional analysis of the system in the sharp interface 
limit A <C 1, as performed in the following section. 

The electrolyte chemical potential fie is an experimen- 
tally accessible independent parameter, as it corresponds 
to the applied potential of the system, forcing Li insertion 
or extraction to occur. Thus, it is important to consider 
the dependence of the wave width and velocity on this 
parameter, as a systematic study of this dependence may 
eventually lead to an experimentally feasible method of 
testing the SRL model. Fig. [6] shows the numerically 
determined dependence of w and v on p,e. As shown in 
Fig. [5^; the wave width exhibits a weak; nonlinear depen- 
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FIG. 5: Scaling of computed wave (a) width and (b) velocity 

with A. Parameters: d = 5, k, — 1, jle — 0.5, c(i, 0) — 
0.1 + 0.8exp(-x2). 



dence on fig, embodied in the scaling function i^^, that 
arises in the dimensional analysis described in the follow- 
ing section. The velocity dependence in Fig. [61d shows 
how the system transitions from extraction waves to in- 
sertion waves as jle increases from negative to positive. 
Large negative values of jle strongly force Li removal, re- 
sulting in extraction waves with large velocities. As fie 
increases, the extraction wave velocity declines steeply, 
until zero velocity is obtained at fie ~ —0.5. This point 
is the transition from extraction to insertion. Insertion 
waves with increasing velocities are produced as fie in- 
creases beyond the transition point. Analogously to the 
width, the velocity dependence is contained within a scal- 
ing function F„ discussed in the following section. Note 
that the width and velocity profiles are not symmetric 
about the transition point, and the minimum width does 
not correspond to zero velocity. The asymmetry in the 
width and velocity result from the asymmetry in the ho- 
mogeneous net reaction rate i?ins — ^oxt- As described 
previously, i?ins — ^^cxt must have three roots over the 



composition range < c < 1 in order for phase sepa- 
ration to occur. This requirement imposes a restricted 
range on fig for a given set of parameters. The extreme 
negative and positive values of fie in Fig. [S] are there- 
fore fixed; there are no traveling wave solutions beyond 
them. The bounding values of fie can be determined an- 
alytically by solving (|43p for the fif corresponding to the 
extrema g± given by (|44[) . as these compositions define 
the limits of the phase separation range. We obtain 

/ 3/2 \ 

A? = 5(1 - 2.g±) + In i J + In (48) 

where /ijr is the minimum allowable potential for extrac- 
tion waves, and is the maximum allowable potential 
for insertion waves. The notational ± signs oi g± and fl^ 
are reversed since g- is the minimum extrcmum at which 
there is almost only insertion, hence corresponding to the 
maximum allowable potential and conversely for g+ 
and yUjT. The limitations on fie are physically meaning- 
ful. In principle, transport can be driven by an arbi- 
trarily positive or negative fie- However, such a strong 
fie effectively increases the overall surface reaction rate k 
and pushes the system out of the SRL transport regime 
Da ^ 1. Indeed, for sufficiently fast surface reactions, 
the system becomes a BTL process {Da ^ 1) with phase 
transformation governed by shrinking core type dynam- 
ics. 

Not all initial conditions give rise to traveling waves, as 
shown in Fig. [71 The key requirement for the formation 
of traveling waves is that the initial condition supports 
both addition and removal of Li in the domain, that is 
i?(x, 0) must change sign. The simultaneous addition and 
removal of material sharpen the initial composition fluc- 
tuation to a phase separating wavefront. In the failed 
insertion event shown in Fig. [7^, only extraction occurs, 
and the initial composition perturbation decays to a uni- 
form concentration of gi. Similarly, Fig. [Tja presents a 
failed extraction event where the initial composition de- 
pression fills up to a uniform concentration of g^. Note 
that R{x, 0) depends on the Laplacian of the initial con- 
centration profile c{x,0), and thus different initial con- 
ditions with equal composition ranges but varying spa- 
tial distributions will or will not produce traveling waves. 
This behavior has been numerically verified. 

The dependence of the formation of traveling waves on 
the initial condition relates to the nucleation of the phase 
separation. During insertion, for example, the lithiated 
phase may first nucleate at some atomic scale inhomo- 
geneity on the crystal surface where it is energetically 
favorable for Li atoms to collect. Our theory does not 
account for such features, and therefore it cannot model 
the initiation of the phase separation. Moreover, as the 
FeP04/LiFeP04 phase boundary width A is nearly at 
the atomic scale, applying a continuum equation for nu- 
cleation below this scale may not be physically relevant. 
We note, however, that continuum nucleation could be 
studied in other SRL systems with larger A. 
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FIG. 6: Dependence of computed wave (a) width and (b) 
velocity on fie. Parameters: 5 = 5, k = 1, A = 1, c(x,0) — 
giifle) + (l/2)(p3(Ae) - gi{iie)){ta.nh{x) + 1). 



E. Scaling with material constants 



FIG. 7: Numerically computed failed insertion and extraction 
events. Dashed lines show initial conditions; solid lines show 
concentration profiles at uniformly spaced times, with arrows 
indicating direction of profile movement. For both insertion 
and extraction: 5 = 5, ft = 1, A = 1. (a) Failed insertion 
event with /ie = 0.5, c(i,0) = 0.1 + O.lexp(-f^). (b) Failed 
extraction event with /le — —1, c(5;,0) = 0.9 — 0.1exp(— i^). 



A complete characterization of wave dynamics in the 
SRL regime is beyond the scope of this paper, but we 
conclude this section by summarizing some key features 
from the analysis above, with dimensions restored. 

Wave nucleation. Starting from a pure phase and ad- 
justing the external chemical potential in the electrolyte 
fie, or the interaction strength a, the nucleation of a 
phase transition wave occurs whenever a random concen- 
tration fluctuation produces a region of sufficiently high, 
stable concentration of the opposite phase. The gradi- 
ent penalty then sharpens the wave, and it propagates 
by the addition or removal of ions at the wavcfront, due 
it its elevated chemical potential compared to both pure 
phases. 



the traveling wave profile is given by 

W — LF. ( — ° Me fccxt \ ^^g^ 

\ L ' pkT ' pkT ' fcins / 

where Fw{X,a, fie, k) is a scaling function. Since the 
width w should not depend on the size of the crys- 
tal L in the limit of a sharp phase boundary, we have 
w = F-uj ~ A/u, or w oc A for A <C 1, consistent with the 
numerical solutions above in Fig.[S^. With units restored, 
wc see that the width is set by the phase-boundary thick- 
ness 



Wave width. By dimensional analysis, the width w of 



as expected, although it may also depend weakly on a, 
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/ie, and K. A slice of the fig dependence is shown in 
Fig. Hi. 

Wave speed. The speed of the traveling wave has a 
similar form. 



^ f ^ ^ Me fccxt 

V = —Fy 



T \ L ' pkT ' pkT ' fcj, 



(51) 



where Fy{X,a, p,e, k) is another scaling function. Once 
again, the limit of a sharp phase boundary requires, v = 
Fy ~ A/« or w oc A/r for A ^ 1, consistent with the 
numerical solutions in Fig. [SJs. 

Recalling the time unit, we obtain a general expression 
for the wave speed in the limit of a sharp phase boundary. 



V pkT pLy 



He 



pkT' pkT' fci, 



(52) 



Note that the wave speed decreases with increasing crys- 
tal thickness, u oc \/Ly since it takes longer for reactions 
to fill each bulk channel, as the SRL phase transforma- 
tion proceeds sweeps across the crystal. The speed is 
also proportional to the gradient penalty coefficient K 
and the insertion rate constant fcins- It also should de- 
crease with the strength of the interaction between ions 
in the crystal a, which drives phase separation. The wave 
velocity can also be controlled externally by varying the 
chemical potential of ions in the electrolyte /Ze, as shown 
in Fig. [5)3. 



IV. RESPONSE TO AN APPLIED VOLTAGE 

A. BTL versus SRL dynamics 

Transport in electrode materials is often studied by 
measuring the current response of the material to an ap- 
plied potential. In an isotropic BTL process, a potential 
step induces a current proportional to the diffusion lim- 
ited flux of ions across the electrode/electrolyte interface. 
The response time of any linear or nonlinear diffusion 
limited process, such as assumed in the shrinking core 
model, is given by the characteristic time tn — /D. 
The flux for small times, and hence the current, can be 
found analytically from the similarity solution for diffu- 
sion in a semi-infinite domain. The resulting expression 
for the current is known as the Cottrell equation. 



'Cottrell 



neAp\ 



(53) 



where n is the number of electrons transferred and A is 
the electrode particle area. The Cottrell current response 
forms the basis of the Potentiostatic Intermittent Titra- 
tion Technique (PITT) that is commonly used to measure 
the diffusivity of materials [131 ■ However, cathodes that 
operate in a SRL transport regime where bulk diffusion 
in a preferred direction is fast relative to ion transfer at 



the electrode/electrolyte interface cannot be assumed to 
follow Cottrell dynamics. 

In our model for single crystal LiFcP04, the chemi- 
cal potential of the electrolyte p,e serves as the applied 
potential to the system. For an appropriate p,e, sharply 
defined waves of Li propagate across the crystal surface 
as Li transfer occurs. A composition fluctuation initiates 
each wave, and all fully developed waves in the system 
travel with the same, constant velocity, in a flat plate-like 
particle. The total flux of ions across the two xz surfaces 
of the crystal is determined by the integral of the net in- 
sertion rate, so the current response of a single crystal is 
given by 



■// 



PsRdx dz, 



(54) 



Note that since the reaction rate R is zero at the equilib- 
rium phase compositions gi and 53, only localized wave- 
fronts spanning these compositions contribute to the cur- 
rent. 

The scaling of the SRL response is radically different 
from that of the Cottrell BTL response. Ignoring ge- 
ometrical effects and nucleation (discussed below), the 
basic scaling of the response time for a single crystal is 
given by 



Lt 



_L _ 
V A 



= L 



pkT pLy 
K 2psku 



(55) 



Note that the response time is proportional to two geo- 
metrical lengths, the depth of the channels Ly and the 
length L over which the waves are propagating. More 
importantly, the time is determined by the surface reac- 
tion rate kj^g and not by the bulk diffusion coefficient D. 
The ratio of the time scales for BTL dynamics and SRL 
dynamics is an effective Peclet number, 



(56) 




pkT pLyD 



which measures the importance of wave propagation at 
the diffusive time scale. However, this is once again the 
Damkohler number, since the reaction time is set by wave 
propagation, tn = ty, in the SRL regime. 



B. Plate-like crystals 

To develop a general picture of SRL phase transfor- 
mation dynamics in a rechargeable battery cathode, we 
first consider the case of flat plate-like crystals of con- 
stant depth Ly analyzed in the previous section. A fully 
developed wavefront moving with constant velocity sup- 
ports a steady current, and there is a sudden loss in the 
current when two wavefronts merge and are replaced by 
an equilibrium composition. Fig. [5] shows this declining 
staircase form for the current in a system with two im- 
pinging waves. The spike in the current at the time of 
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FIG. 8; Numerically computed phase transformation dynam- 
ics and flux response of two impinging insertion waves, (a) 
Concentration profiles of impinging waves. Dashed line shows 
initial condition; solid lines show concentration profile at var- 
ious times during impingement, (b) Flux response of imping- 
ing waves. Parameters: a = 5, k = 1, A = 1, /ie = 0.5, 
c{x, 0) = (71 + (33 - 5i)[exp(-(i -f 7)'') + exp(-(i - 7)^)]. 



collision is due to the gradient penalty term in the reac- 
tion rate acting on the sharp composition profile of the 
merging waves. We note that the magnitude of the spike 
is large in this simulation since there are only two waves; 
in an actual system with many waves, the surge in cur- 
rent from any individual collision would be small relative 
to the total current being sustained. 

Thus, the current response of LiFeP04 is governed by 
the overall rate of its transformation through concurrent 
nucleation and growth of waves. Phase nucleation likely 
occurs by both heterogeneous and homogeneous mecha- 
nisms. Recent first principles computations have found 
that the chemical potential of Li varies considerably over 
the surface of the equilibrium crystal shape [2^. Con- 
sequently, different crystal faces may be energetically fa- 
vorable for heterogeneous phase nucleation during Li in- 



sertion and extraction. 

An example of a single crystal undergoing heteroge- 
neous and homogeneous nucleation and growth is illus- 
trated in Fig. O Fig. [^^ shows xy cross sections of the 
crystal at a sequence of successive times ti,...^tQ, and 
Fig.Hb presents the corresponding profile of the total cur- 
rent /. At time ti, heterogeneous nucleation at the crys- 
tal edges has produced two fully developed wavefronts, 
each moving with constant velocity v and sustaining a 
normalized current of unity. Therefore, the crystal sup- 
ports the total current / = 2 at this time. At some 
time between ti and t2, heterogeneous nucleation occurs 
at the two rightmost surface defects of the crystal, in- 
dicated by notches. Once these nucleation events grow 
into fully developed waves, there are 6 propagating wave- 
fronts carrying a total current of / = 6. The rightmost 
waves merge at time such that 4 wavefronts are de- 
stroyed, and consequently, only two traveling wavefronts 
remain and the current drops to / = 2. Homogeneous nu- 
cleation at some location in the untransformed fraction 
of the material occurs at time ^4, and the two additional 
wavefronts created increase the current to / = 4. The 
rightmost waves combine at time is, and as most of the 
material is transformed and only two wavefronts remain, 
the current again drops to / = 2. Finally, at time t^, the 
material is fully transformed and can no longer sustain a 
current. 



C. Other crystal shapes 

The wave dynamics can depend sensitively on the crys- 
tal shape in the SRL regime. This general fact can be eas- 
ily seen from our analysis of a plate-like crystal: the wave 
velocity ([5^ depends on the local depth Ly{x,z) of the 
fast diffusion channels in the bulk, as well as the local sur- 
face orientation, through the surface-site density Ps{x, z). 
The analysis of SRL phase-transformation dynamics for 
arbitrary crystal shapes is a challenging problem, left for 
future work. 

Here, we simply indicate how scalings by considering 
the limit of slowly varying depth, Ly(x), and assuming 
the ID wave dynamics for a fiat surface arc only slightly 
perturbed. The local wave velocity is then 



(57) 



where = 2ps{xw)kins\/K/ [p^kT] sa constant. This 
ordinary differential equation can be solved for the po- 
sition of the wave Xw{t) from a given wave nucleation 
event for any slowly varying shape Lyixw)- For exam- 
ple, consider the case of a cylinder, Ly{x) 



(ignoring that the shape is not slowly varying at the 
ends). In that case, the equation can be solved analyti- 
cally, e.g., for nucleation at one edge x = —R. For early 
times, the wave velocity decays from its initial value as 
dxw/dt oc Dti^t"^/^ , due to the increasing depth of the 
bulk channels. 
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crystal thus scales with the total volume (in contrast to 
diffusive BTL dynamics, where the time scales like the 
cross-sectional area). Physically, ions are being inserted 
or extracted at roughly a constant rate, since the phase 
boundary is assumed to have a constant exposed length 
at the surface for ID dynamics. 

For spheres and other 3D shapes, the wavefront will not 
remain flat, and the full 2D depth-integrated dynamics 
will need to be solved with variable Ly{x, z) and Ps{x^ z). 
However, the current may remains roughly constant dur- 
ing wave propagation, since the filling time for a channel 
is proportional to its length, at constant surface reaction 
rate. In that case, the results of the previous section 
for current versus time in flat plate-like single crystals 
may not be substantially modified with more complicated 
shapes, although statistical fiuctuations due to random 
nuclcation events will be different. 



^^^^ 



D. Composite cathode response 




FIG. 9: Schematic diagram of overall phase transformation 
and current response of a single plate-like crystal undergoing 
SRL transport, (a) Sequence of xy cross sections of crystal at 
times t\,. . . ,1^, illustrating phase transformation of material 
through concurrent nucleation and growth of traveling waves. 
Each fully developed wavefront moves with velocity v and 
sustains a constant current of unity, (b) Declining staircase 
current response of crystal, with labeled times corresponding 
to illustrations in (a). 



In spite of variable wave speed, however, the current 
remains constant during wave propagation in this ap- 
proximation, as in the case of a flat plate, regardless of 
the crystal shape. The reason is that wave propagation 
at speed v{xu]) engulfs channels of length Ly(xw), so the 
total current, proportional to vLy, remains constant ac- 
cording to (j57p . The time for a wave to engulf the entire 



It is important to note that Fig. [9] represents only one 
possible realization of the transformation and current re- 
sponse of the crystal. Heterogeneous nucleation may oc- 
cur at different edges or surface defects at different times 
for different insertion and extraction cycles. Homoge- 
neous nucleation would be spatially distributed in some 
random fashion. Therefore, to determine the overall cur- 
rent response of a composite cathode composed of many 
individual crystals, we must consider the statistical dis- 
tributions of the nucleation events. 

For homogeneous nucleation, we may assume that 
the nucleation rate is uniform across the crystal sur- 
face. Nucleation events in the untransformed material 
are independent, and the presence of a previously nucle- 
ated wave does not influence the likelihood of nucleation 
around that wave. With these assumptions, the nucle- 
ation events arc distributed as a Poisson process in time, 
and we may invoke the Johnson-Mchl-Avrami equation 
for the overall transformation rate of the material 



-Gt" 



(58) 



where ^ is the untransformed fraction of the material, G 
is a factor dependent on the dimensionality of the pro- 
cess, and the integer n > 1; for a ID line nucleation 
process, Gt" ~ t^. The Johnson-Mehl-Avrami equation 
therefore specifies that the overall transformation rate of 
the material follows a sigmoidal shape. We thus expect 
that the current also follows this sigmoidal response for 
concurrent homogeneous nuclcation and growth. 

In the case of heterogeneous nucleation, we may as- 
sume that for many defects in many particles, the het- 
erogeneous nucleation events are distributed as a Poisson 
process in space. Thus, the qualitative form of the av- 
erage response is expected to follow the same sigmoidal 
shape as homogeneous nucleation, assuming all electrode 
particles are at the same potential (driving force), which 
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FIG. 10: Schematic illustration of composite cathode re- 
sponses for Cottrell (BTL) and SRL phase transformation 
mechanisms. 

may only be true under low rate conditions or for thin 
electrodes [llli]. 

Thus, we have found that individual crystals have a 
declining staircase form for their current response to an 
applied electrolyte potential fj,e if a substantial number 
of nuclcation events occur on the timescale of complete 
transformation. For many crystals, we may consider that 
the homogeneous and heterogeneous nucleation events 
are Poisson processes in time and space, respectively. 
The overall transformation and current response of the 
material is then given by the average of many distinct 
staircase responses, resulting in a characteristic sigmoidal 
curve 

J-e-<^*", (59) 

that is strikingly different than the Cottrell response that 
is commonly assumed. Fig. [TU] compares the sigmoidal 
and Cottrell responses. We conclude that PITT mea- 
surements in material undergoing SRL transport do not 
measure the diffusivity. Rather, these measurements pro- 
vide some measure of the kinetic parameters in the sur- 
face reaction rates that are controlling the overall rate of 
the material transformation. 



E. Discussion 

The model we propose is based on a few key experi- 
mental and theoretical findings. (1) Li migrates rapidly 
in the (010) direction {y in Fig. [1]) creating either filled 
or empty channels that completely penetrate the mate- 
rial. This makes it possible to coarse-grain to a two- 
dimensional model in which the surface concentration 
defines the concentration in a channel. (2) A linear inter- 
face exists on the particle surface between filled and un- 
filled sites, and growth of one phase at the expense of the 



other occurs by displacement of this interface. Such one- 
dimensional growth is supported by experimental obser- 
vations [l^l and by a recent Johnson-Mehl Avrami analy- 
sis of the growth exponents [ssj . Since one would expect 
two-dimensional growth on the surface for an isotropic 
material, the one-dimensional growth has to find its ori- 
gin in the crystallography of the material. Undoubtedly, 
it is either the anisotropy of the interfacial strain en- 
ergy - because of different coherence strains or elastic 
constants or the anisotropy of the interfacial energy 
which causes the system to prefer a single interface plane. 

In our model, the transformation rate is controlled by 
transfer of the ions from the electrolyte to the (010) sur- 
face. Since during growth insertion only takes place at 
the interface, the energy of the Li ions at this interface is 
a key quantity in determining the effective transfer rate 
from the solution. Unlike in a core-shell model where 
the ability of the particle to take up current declines as 
(de)lithiation proceeds, in our model the transformation 
rate is nominally constant, until the particle is either fully 
transformed, new nucleation events occur, or two wave 
fronts impinge. Such transformation kinetics is funda- 
mentally different from Cottrell-like behavior. 

How this transformation kinetics manifests itself in 
the observable voltage-current response may depend very 
much on the structure of the macroscopic electrode in 
which the active LiFeP04 particles are embedded. In ad- 
dition to the active material, a typical electrode contains 
about 5-10 wt% polymeric binder and 5-15 wt% carbon 
black to enhance electronic conductivity through the elec- 
trode. Some porosity is also created in the electrode to 
allow the electrolyte to penetrate and transport the Li+ 
ions to and from the active material. If the conductive 
pathways for the Li+ ions and electrons are sufficient, 
all particles will be at the same potential and experience 
a similar driving force for transformation. Under these 
conditions, and assuming stochastic nucleation, we ex- 
pect that the overall current response to a potential pulse 
is sigmoidal. Recent work indicates that such equipoten- 
tial conditions across the electrode only apply at rather 
low charge and discharge rates, or for very thin electrodes 
[3^ . If such electrical resistance along the thickness of 
the electrode plays a role, the collective current response 
of the system could be viewed as of a sum of sigmoidals, 
each with a different driving force, but time-dependent 
screening effects would also need to be taken into ac- 
count. 



V. CONCLUSION 

We have proposed a general continuum theory for Li 
transport in single crystal LiFeP04 based on a thermo- 
dynamically sound phase field formulation of the free en- 
ergy. In various limits of the characteristic timescales 
for bulk and surface transport, this theory captures the 
shrinking core and other BTL models. In the limit where 
fast bulk diffusion in ID equilibrates the bulk concen- 
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tration to the surface concentration, our model gives a 
new type of equation describing SRL phase transition 
dynamics. We find that this model exhibits traveling 
wave solutions that qualitatively agree with the experi- 
mental observations. The main implication of our work 
for battery modeling is that the current response of SRL 
systems, such as LiFeP04, are not governed by the clas- 
sical, bulk diffusion limited Cottrell model that is com- 
monly assumed. Our work also focuses attention to the 
importance of the Li"*" and electron delivery to the proper 
surface of LiFeP04 in order to achieve fast charge absorp- 
tion. While much effort in the experimental literature 
has focused on electron delivery (e.g. by carbon coating 
or conductive Fe2P contributions) [3J,l35|, little empha- 



sis seems to have been place on rapid transport of Li+ 
towards the surface where it can penetrate. Finally, we 
note that the SRL model developed here may be appli- 
cable in other materials where surface transfer effects are 
rate limiting, such as nanoporous materials. 
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